# -*- coding: utf-8 -*-
"""
Created on Sat Jan 26 08:50:06 2019

@author: Salem
"""

import numpy as np
import numpy.linalg as la
import numpy.random as npr


import scipy.optimize as op

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from stl import mesh

#-------------------------------------
from pathlib import WindowsPath
from pathlib import Path
import os

cwd = Path(os.getcwd())

Dropbox = Path(__file__).parents[3] 

mesh_dir = Dropbox / "With Arkhat"  / "Meshes" 
results_dir = cwd  / "Run Results"
top_dir = mesh_dir / "Tops A" 
top_open_dir = mesh_dir / "Tops Open" 
top_CH_dir = mesh_dir / "Tops CH" 

tip_open_dir = mesh_dir / "Tips Open" 

wtr_dir = mesh_dir / "Watertight B"

full_dir = mesh_dir / "Full Skulls STL" 
cuts_dirA = mesh_dir / "Cut Beaks"
cuts_dirB = mesh_dir / "Cut Beaks Aligned"

other_full_dir = mesh_dir / "Other Species" 
#-------------------------------------


#names of beaks
beak1 =  "G.FuliginosaA"
beak2 =  "G.FuliginosaB"
beak3 =  "G.FuliginosaC"
beak4 =  "G.FuliginosaD"
beak5 =  "G.FuliginosaE"
beak6 = "G.FuliginosaF"
beak7 =  "G.FortisA"
beak8 =  "G.FortisB"
beak9 =  "G.FortisC"
beak10 =  "G.FortisD"
beak11 =  "G.FortisE"
beak12 = "G.MagnirostrisA"
beak13 = "G.MagnirostrisB"
beak14 = "G.MagnirostrisC"
beak15 = "G.MagnirostrisD"
beak16 = "G.MagnirostrisE"
beak17 =  "G.ConirostrisA"
beak18 =  "G.ConirostrisB"
beak19 =  "G.ConirostrisC"
beak20 =  "G.ConirostrisD"
beak21 =  "G.ConirostrisE"
beak22 =  "G.ConirostrisF"
beak23 =  "G.ConirostrisG"
beak24 =  "G.ScandensA"
beak25 =  "G.ScandensB"
beak26 =  "G.SeptentrionalistA"
beak27 =  "G.SeptentrionalistB"
beak28 =  "G.DifficilisA"
beak29 =  "G.DifficilisB"
beak30 =  "C.PallidusA"
beak31 =  "C.PallidusB"
beak32 = "C.PallidusC"
beak33 = "C.PallidusD"
beak34 =  "C.ParvulusA"
beak35 =  "C.ParvulusB"
beak36 =  "C.ParvulusC"
beak37 =  "C.ParvulusD"
beak38 =  "C.ParvulusE"
beak39 =  "C.PsittaculaA"
beak40 = "C.PsittaculaB"
beak41 =  "C.PsittaculaC"
beak42 =  "C.PsittaculaD"
beak43 =  "P2.InornataA"
beak44 =  "P2.InornataB"
beak45 = "C2.OlivaceaA"
beak46 =  "C2.OlivaceaB"
beak47 =  "C2.OlivaceaC"
beak48 =  "C2.FuscaA"
beak49 =  "C2.FuscaB"
beak50 =  "P.CrassirostrisA"
beak51 =  "P.CrassirostrisB"
beak52 =  "P.CrassirostrisC"
beak53 =  "P.CrassirostrisD"
beak54 =  "P.CrassirostrisE"
beak55 =  "T.BicolorA"
beak56 =  "T.BicolorB"
beak57 =  "T.BicolorC"
beak58 =  "T.BicolorD"
beak59 =  "T.BicolorE"


names = np.array([beak1, beak2, beak3, beak4, beak5, beak6, beak7, beak8, beak9, 
                  beak10, beak11, beak12, beak13, beak14, beak15, beak16, beak17, beak18,
                  beak19, beak20, beak21, beak22, beak23, beak24, beak25, beak26, beak27,
                  beak28, beak29, beak30, beak31, beak32, beak33, beak34, beak35, beak36, beak37,
                  beak38, beak39, beak40, beak41, beak42, beak43, beak44, beak45, beak46, beak47,
                  beak48, beak49, beak50, beak51, beak52, beak53, beak54, beak55, beak56, beak57, beak58, beak59]) 


other_beak1 = "A.Lathami"
other_beak2 = "B.ruficapillus"
other_beak3 = "E._helias"
other_beak4 = "X.zosterops"

other_names =  np.array([other_beak1, other_beak2, other_beak3, other_beak4]) 


#=========================================================================================================================
#Takes and stl file and gets the connectivity data: vertex indices and where they appear in faces and edges
#========================================================================================================================= 
def get_connectivity(mesh_name):
    
    stl_mesh = mesh.Mesh.from_file(mesh_name)
    
    points  = stl_mesh.points.round(decimals=4)
    
    vertices, inverts  = np.unique(points.reshape((points.size//3, 3)), axis=0, return_inverse=True)
    
    faces = inverts.reshape((inverts.size//3, 3)) #indices of triangles. 
    
    #edges = np.unique(faces[:, [[0,1], [0,2], [1,2]]].reshape((faces.size, 2)), axis=0)
    
    
    return  vertices, faces
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#========================================================================================================================= 



#=========================================================================================================================
# takes the vertices and faces and generates then saves an stl mesh
#========================================================================================================================= 
def make_stl(vertices, faces, mesh_location):
    
    # Create the mesh
    stl_mesh = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
    for i, f in enumerate(faces):
        for j in range(3):
            stl_mesh.vectors[i][j] = vertices[f[j],:]


    print("saving to: ", mesh_location, "\n\n")

    # Write the mesh to file mesh_name
    stl_mesh.save(mesh_location)
    
    return 
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#========================================================================================================================= 


#=========================================================================================================================
    # chooses a point on a sphere embedded in ndim (dimensions of sphere is ndim-1)
#=========================================================================================================================      
def sample_spherical(npoints, ndim=3):
        vec = np.random.randn(ndim, npoints)
        vec /= np.linalg.norm(vec, axis=0)
        return vec.transpose()
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================   
    

#=========================================================================================================================
# gets the moment of inertia assuming unit mass (could be changed easily)
#========================================================================================================================= 
def get_moment_tensor(vertices):
    
    return np.einsum('ai, ai', vertices, vertices)*np.eye(3) - np.einsum('ai, aj -> ij', vertices, vertices)
   
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================    
    
    
#=========================================================================================================================
# returns the eigendirections of the moment tensor. each direction is a row in the matrix
#========================================================================================================================= 
def get_moment_directions(vertices):
    
    moment_tensor = get_moment_tensor(vertices)
    
    return np.transpose(la.eigh(moment_tensor)[1]) #gives the principle directions
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================    
    
#=========================================================================================================================
# align principle axes with the xyz
#========================================================================================================================= 
def align_moment_direction(vertices):
    
    orthogonal_mat = get_moment_directions(vertices).transpose()
    
    #multiply an extra direction so that an aligned mesh stays the same
    reflection = np.diag(np.sign(np.diagonal(orthogonal_mat)))
    
    orthogonal_mat = np.dot(orthogonal_mat, reflection)
    
    plane_verts = np.array([[0, -1, -1], [0, -1, 1], [0, 1, 1],[0, 1, -1]])
    plane_faces = np.array([[0, 1, 3], [1, 2, 3]])
    
    #gives the rotated vertices 
    return  np.dot(vertices, orthogonal_mat), (np.dot(plane_verts, orthogonal_mat), plane_faces)
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================     


#=========================================================================================================================
# gets the plane aligned with the axes of the mesh, will be used for cutting along appropriate axes
#========================================================================================================================= 
def get_cutting_plane(mesh_name,  plane_name, isSaveSTL=True):
    
    vertices, faces = get_connectivity(mesh_name)
    
    #no transpose in this case since the plane should transform the opposite way to the mesh when aligned
    orthogonal_mat = get_moment_directions(vertices)
    
    
    #make the rotated plane
    plane_verts = np.array([[0, -1, -1], [0, -1, 1], [0, 1, 1],[0, 1, -1]])
    plane_verts = np.dot(plane_verts, orthogonal_mat)
    plane_faces = np.array([[0, 1, 3], [1, 2, 3]])
    
    if isSaveSTL:
        make_stl(plane_verts, plane_faces, plane_name)
        return
    #gives the rotated vertices 
    return  plane_verts, plane_faces
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================     



#=========================================================================================================================
# takes in an stl mesh, alings the principal axes of inertia to xyz, then overwrites the new mesh to the file.
#========================================================================================================================= 
from scipy.spatial import ConvexHull
def make_aligned_stl(mesh_name, load_dir=full_dir, save_dir=full_dir, origin_to_geometry=True):
    
    vertices, faces = get_connectivity(load_dir / (mesh_name + '.stl'))
    
    if origin_to_geometry:
        new_verts = vertices - vertices.mean(axis=0)
    
    new_verts, ignore = align_moment_direction(new_verts)
    
    make_stl(new_verts, faces, save_dir / (mesh_name + '.stl'))
    
    '''
    if isPlot: plot(new_verts)
    
    if isSaveHull:
        
        hull = ConvexHull(new_verts)
        
        hull_verts = hull.points
        
        #shift vertices so that the max_position = - min_position
        hull_verts[:,0] -= 0.5*(hull_verts[:,0].max() + hull_verts[:,0].min()) # x-direction
        hull_verts[:,1] -= 0.5*(hull_verts[:,1].max() + hull_verts[:,1].min()) # y-direction
        hull_verts[:,2] -= 0.5*(hull_verts[:,2].max() + hull_verts[:,2].min()) # z-direction
        
        print(hull_verts[hull.vertices].shape)
        
        plot(hull_verts[hull.vertices])
        
        make_stl(hull_verts, hull.simplices, hull_name)    
    '''    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
   
#=========================================================================================================================
# takes the unaligned meshes and alings them for all sample names
#========================================================================================================================= 
def align_all_stls(mesh_names = names, load_dir=cuts_dirA, save_dir=cuts_dirB):
    
    for name in mesh_names:
        make_aligned_stl(name, load_dir=load_dir, save_dir=save_dir)
        
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================     

#=========================================================================================================================
# takes in an stl mesh and saves the convex hull after centering
#========================================================================================================================= 
def make_centered_hull(mesh_name,  hull_name=''):
    
    vertices, faces = get_connectivity(mesh_name) 
        
    hull = ConvexHull(vertices)
        
    hull_verts = hull.points
        
        #shift vertices so that the max_position = - min_position
    hull_verts[:,0] -= 0.5*(hull_verts[:,0].max() + hull_verts[:,0].min()) # x-direction
    hull_verts[:,1] -= 0.5*(hull_verts[:,1].max() + hull_verts[:,1].min()) # y-direction
    hull_verts[:,2] -= 0.5*(hull_verts[:,2].max() + hull_verts[:,2].min()) # z-direction
        
    print(hull_verts[hull.vertices].shape)
        
    plot(hull_verts[hull.vertices])
        
    make_stl(hull_verts, hull.simplices, hull_name)    
        
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
    
#=========================================================================================================================
# takes in an stl mesh and saves the convex hull after centering
#========================================================================================================================= 
def make_hull(mesh_name, load_dir=top_dir, save_dir=top_CH_dir):
    
    vertices, faces = get_connectivity(load_dir / (mesh_name + '.stl')) 
        
    hull = ConvexHull(vertices)
        
    hull_verts = hull.points
        
   
    print(hull_verts[hull.vertices].shape)
        
    #plot(hull_verts[hull.vertices])
        
    make_stl(hull_verts, hull.simplices, save_dir / (mesh_name + '.stl'))    
        
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
    
#=========================================================================================================================
# takes in an stl mesh, alings the principal axes of inertia to xyz, then overwrites the new mesh to the file.
#========================================================================================================================= 
def make_centered_mesh(mesh_name):
    
    vertices, faces = get_connectivity(mesh_name) 
    
    
        
        #shift vertices so that the max_position = - min_position
    vertices[:,0] -= 0.5*(vertices[:,0].max() + vertices[:,0].min()) # x-direction
    vertices[:,1] -= 0.5*(vertices[:,1].max() + vertices[:,1].min()) # y-direction
    vertices[:,2] -= 0.5*(vertices[:,2].max() + vertices[:,2].min()) # z-direction
        
        
    plot(vertices)
        
    make_stl(vertices, faces, mesh_name)    
        
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
    

#=========================================================================================================================
# takes in an stl mesh, then shifts the tip to the origin
#========================================================================================================================= 
def shift_tip_to_origin(mesh_name, load_dir=top_open_dir, save_dir=top_open_dir):
    
    vertices, faces = get_connectivity(load_dir / (mesh_name + '.stl')) 
    
    X = vertices[:,0]
    
    X_tip  = np.amin(X)
    
    index_of_tip = np.where(X==X_tip)[0][0]

    position_of_tip = vertices[index_of_tip] 
    
    #shift vertices so that the min_position = 0
    vertices -= position_of_tip
    
   #plot(vertices)
        
    make_stl(vertices, faces, save_dir / (mesh_name + '.stl'))    
        
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
    

#=========================================================================================================================
# apply a transformation to a mesh for comparison
#========================================================================================================================= 
def apply_transformation(mesh_name, transformation, translation):
    
    vertices, faces = get_connectivity(mesh_name)
    
    #center_of_mass = np.mean(vertices, axis=0)
    new_verts = vertices #- center_of_mass
    
    #gives the transformed vertices 
    new_verts = np.dot(new_verts, transformation) + translation
       
    #shift vertices so that the min_position = 0
    #new_verts[:,0] -= new_verts[:,0].min() # x-direction
    
    make_stl(new_verts, faces, mesh_name)
    
    return  
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================     

    

#=========================================================================================================================
    # plot the vertices as a scatter plot
#========================================================================================================================= 
def plot(vertices):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    x =vertices[:, 0]
    y =vertices[:, 1]
    z =vertices[:, 2]



    ax.scatter(x, y, z, c='r', marker='o')

    ax.set_xlabel('X Label')
    ax.set_ylabel('Y Label')
    ax.set_zlabel('Z Label')
    
    return
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================      
    


    
#=========================================================================================================================
# finds a list of lengths, widths and heights for a set of meshes
#=========================================================================================================================  
def get_all_dimensions(mesh_names=names, load_dir=tip_open_dir, save_name="dimensions_tip.csv"):
    
    num_meshes = mesh_names.shape[0]
    
    dimensions = np.zeros((num_meshes, 3))
    
    for i, mesh_name in enumerate(mesh_names):

    
        verts, faces  = get_connectivity(load_dir /(mesh_name + '.stl'))
       
        
        dimensions[i] = get_dimensions(verts)
    
    plot(dimensions)
    
    #norms = (la.norm(dimensions, axis=1).reshape((num_meshes, 1)))
    #plot(dimensions/norms)
    
    np.savetxt(results_dir / save_name, dimensions, delimiter=',')
    
    return dimensions
    
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  


    

#=========================================================================================================================
# find the length and width and height for a set of vertices
#=========================================================================================================================  
def get_dimensions(vertices):
    
    length = vertices[:,0].max() -  vertices[:,0].min()
    
    width = vertices[:,1].max() -  vertices[:,1].min()
    
    height = vertices[:,2].max() -  vertices[:,2].min()
    
    return [length, width, height]
    
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================  
    



    
skull_name1 =  mesh_dir / "Full Skulls STL" / "G.FortisA.stl"
beak_name1 =  mesh_dir / "Cut Beaks B" / "G.FortisA_Cut.stl"  
CH_name1 =  mesh_dir / "Convex Hulls A" / "G.FortisA_CH.stl"   
WTR_name1 =  mesh_dir / "Watertight B" / "G.FortisA_WTR.stl"   
top_name1 =  mesh_dir / "Tops A" / "G.FortisA.stl" 
fit_name1 =  mesh_dir / "Fits A" / "G.FortisA.stl" 

skull_name2 =  mesh_dir / "Full Skulls STL" / "G.MagnirostrisA.stl"
beak_name2 =  mesh_dir / "Cut Beaks B" / "G.MagnirostrisA_Cut.stl"
CH_name2 =  mesh_dir / "Convex Hulls A" / "G.MagnirostrisA_CH.stl" 
WTR_name2 =  mesh_dir / "Watertight B" / "G.MagnirostrisA_WTR.stl" 
top_name2 =  mesh_dir / "Tops A" / "G.MagnirostrisA.stl" 
fit_name2 =  mesh_dir / "Fits A" / "G.MagnirostrisA.stl" 

skull_name3 =  mesh_dir / "Full Skulls STL" / "C.ParvulusB.stl"
beak_name3 =  mesh_dir / "Cut Beaks B" / "C.ParvulusB_Cut.stl"
CH_name3 =  mesh_dir / "Convex Hulls A" / "C.ParvulusB_CH.stl"   
WTR_name3 =  mesh_dir / "Watertight B" / "C.ParvulusB_WTR.stl" 
top_name3 =  mesh_dir / "Tops A" / "C.ParvulusB.stl" 
fit_name3 =  mesh_dir / "Fits A" / "C.ParvulusB.stl" 


skull_name4 =  mesh_dir / "Full Skulls STL" / "C.PsittaculaA.stl"
beak_name4 =  mesh_dir / "Cut Beaks B" / "C.PsittaculaA_Cut.stl"
CH_name4 =  mesh_dir / "Convex Hulls A" / "C.PsittaculaA_CH.stl"   
WTR_name4 =  mesh_dir / "Watertight B" / "C.PsittaculaA_WTR.stl" 
top_name4 =  mesh_dir / "Tops A" / "C.PsittaculaA.stl" 
fit_name4 =  mesh_dir / "Fits A" / "C.PsittaculaA.stl" 
        
skull_name5 =  mesh_dir / "Full Skulls STL" / "G.FortisB.stl"
beak_name5 =  mesh_dir / "Cut Beaks B" / "G.FortisB_Cut.stl"    
CH_name5 =  mesh_dir / "Convex Hulls A" / "G.FortisB_CH.stl"   
WTR_name5 =  mesh_dir / "Watertight B" / "G.FortisB_WTR.stl"  
top_name5 =  mesh_dir / "Tops A" / "G.FortisB.stl" 
fit_name5 =  mesh_dir / "Fits A" / "G.FortisB.stl" 
    
skull_name6 =  mesh_dir / "Full Skulls STL" / "G.FuliginosaB.stl"
beak_name6 =  mesh_dir / "Cut Beaks B" / "G.FuliginosaB_Cut.stl"    
CH_name6 =  mesh_dir / "Convex Hulls A" / "G.FuliginosaB_CH.stl"
WTR_name6 =  mesh_dir / "Watertight B" / "G.FuliginosaB_WTR.stl" 
top_name6 =  mesh_dir / "Tops A" / "G.FuliginosaB.stl" 
fit_name6 =  mesh_dir / "Fits A" / "G.FuliginosaB.stl" 

skull_name7 =  mesh_dir / "Full Skulls STL" / "C.PallidusA.stl"
beak_name7 =  mesh_dir / "Cut Beaks B" / "C.PallidusA_Cut.stl"    
CH_name7 =  mesh_dir / "Convex Hulls A" / "C.PallidusA_CH.stl"
WTR_name7 =  mesh_dir / "Watertight B" / "C.PallidusA_WTR.stl" 
top_name7 =  mesh_dir / "Tops A" / "C.PallidusA.stl" 
fit_name7 =  mesh_dir / "Fits A" / "C.PallidusA.stl" 

skull_name8 =  mesh_dir / "Full Skulls STL" / "C.PallidusB.stl"
beak_name8 =  mesh_dir / "Cut Beaks B" / "C.PallidusB_Cut.stl" 
CH_name8 =  mesh_dir / "Convex Hulls A" / "C.PallidusB_CH.stl"  
WTR_name8 =  mesh_dir / "Watertight B" / "C.PallidusB_WTR.stl" 
top_name8 =  mesh_dir / "Tops A" / "C.PallidusB.stl" 
fit_name8 =  mesh_dir / "Fits A" / "C.PallidusB.stl" 

skull_name9 =  mesh_dir / "Full Skulls STL" / "C.PallidusC.stl"
beak_name9 =  mesh_dir / "Cut Beaks B" / "C.PallidusC_Cut.stl"    
CH_name9 =  mesh_dir / "Convex Hulls A" / "C.PallidusC_CH.stl"   
WTR_name9 =  mesh_dir / "Watertight B" / "C.PallidusC_WTR.stl" 
top_name9 =  mesh_dir / "Tops A" / "C.PallidusC.stl" 
fit_name9 =  mesh_dir / "Fits A" / "C.PallidusC.stl" 

skull_name10 =  mesh_dir / "Full Skulls STL" / "C.PsittaculaB.stl"
beak_name10 =  mesh_dir / "Cut Beaks B" / "C.PsittaculaB_Cut.stl"
CH_name10 =  mesh_dir / "Convex Hulls A" / "C.PsittaculaB_CH.stl"   
WTR_name10 =  mesh_dir / "Watertight B" / "C.PsittaculaB_WTR.stl" 
top_name10 =  mesh_dir / "Tops A" / "C.PsittaculaB.stl" 
fit_name10 =  mesh_dir / "Fits A" / "C.PsittaculaB.stl" 

skull_name11 =  mesh_dir / "Full Skulls STL" / "C.PsittaculaC.stl"
beak_name11 =  mesh_dir / "Cut Beaks B" / "C.PsittaculaC_Cut.stl"
CH_name11 =  mesh_dir / "Convex Hulls A" / "C.PsittaculaC_CH.stl"   
WTR_name11 =  mesh_dir / "Watertight B" / "C.PsittaculaC_WTR.stl" 
top_name11 =  mesh_dir / "Tops A" / "C.PsittaculaC.stl" 
fit_name11 =  mesh_dir / "Fits A" / "C.PsittaculaC.stl" 

skull_name12 =  mesh_dir / "Full Skulls STL" / "G.FortisC.stl"
beak_name12 =  mesh_dir / "Cut Beaks B" / "G.FortisC_Cut.stl"
CH_name12 =  mesh_dir / "Convex Hulls A" / "G.FortisC_CH.stl"  
WTR_name12 =  mesh_dir / "Watertight B" / "G.FortisC_WTR.stl"  
top_name12 =  mesh_dir / "Tops A" / "G.FortisC.stl" 
fit_name12 =  mesh_dir / "Fits A" / "G.FortisC.stl" 

skull_name13 =  mesh_dir / "Full Skulls STL" / "G.FuliginosaA.stl"
beak_name13 =  mesh_dir / "Cut Beaks B" / "G.FuliginosaA_Cut.stl" 
CH_name13 =  mesh_dir / "Convex Hulls A" / "G.FuliginosaA_CH.stl"  
WTR_name13 =  mesh_dir / "Watertight B" / "G.FuliginosaA_WTR.stl"  
top_name13 =  mesh_dir / "Tops A" / "G.FuliginosaA.stl" 
fit_name13 =  mesh_dir / "Fits A" / "G.FuliginosaA.stl"   

skull_name14 =  mesh_dir / "Full Skulls STL" / "G.FuliginosaC.stl"
beak_name14 =  mesh_dir / "Cut Beaks B" / "G.FuliginosaC_Cut.stl" 
CH_name14 =  mesh_dir / "Convex Hulls A" / "G.FuliginosaC_CH.stl" 
WTR_name14 =  mesh_dir / "Watertight B" / "G.FuliginosaC_WTR.stl"   
top_name14 =  mesh_dir / "Tops A" / "G.FuliginosaC.stl" 
fit_name14 =  mesh_dir / "Fits A" / "G.FuliginosaC.stl"    

skull_name15 =  mesh_dir / "Full Skulls STL" / "G.MagnirostrisB.stl"
beak_name15 =  mesh_dir / "Cut Beaks B" / "G.MagnirostrisB_Cut.stl"  
CH_name15 =  mesh_dir / "Convex Hulls A" / "G.MagnirostrisB_CH.stl" 
WTR_name15 =  mesh_dir / "Watertight B" / "G.MagnirostrisB_WTR.stl"   
top_name15 =  mesh_dir / "Tops A" / "G.MagnirostrisB.stl" 
fit_name15 =  mesh_dir / "Fits A" / "G.MagnirostrisB.stl" 

skull_name16 =  mesh_dir / "Full Skulls STL" / "G.MagnirostrisC.stl"
beak_name16 =  mesh_dir / "Cut Beaks B" / "G.MagnirostrisC_Cut.stl" 
CH_name16 =  mesh_dir / "Convex Hulls A" / "G.MagnirostrisC_CH.stl"   
WTR_name16 =  mesh_dir / "Watertight B" / "G.MagnirostrisC_WTR.stl"    
top_name16 =  mesh_dir / "Tops A" / "G.MagnirostrisC.stl" 
fit_name16 =  mesh_dir / "Fits A" / "G.MagnirostrisC.stl" 
  
skull_name17 =  mesh_dir / "Full Skulls STL" / "G.ScandensB.stl"
beak_name17 =  mesh_dir / "Cut Beaks B" / "G.ScandensB_Cut.stl"    
CH_name17 =  mesh_dir / "Convex Hulls A" / "G.ScandensB_CH.stl"  
WTR_name17 =  mesh_dir / "Watertight B" / "G.ScandensB_WTR.stl"    
top_name17 =  mesh_dir / "Tops A" / "G.ScandensB.stl" 
fit_name17 =  mesh_dir / "Fits A" / "G.ScandensB.stl" 

skull_name18 =  mesh_dir / "Full Skulls STL" / "G.ScandensA.stl"
beak_name18 =  mesh_dir / "Cut Beaks B" / "G.ScandensA_Cut.stl"    
CH_name18 =  mesh_dir / "Convex Hulls A" / "G.ScandensA_CH.stl"
WTR_name18 =  mesh_dir / "Watertight B" / "G.ScandensA_WTR.stl" 
top_name18 =  mesh_dir / "Tops A" / "G.ScandensA.stl" 
fit_name18 =  mesh_dir / "Fits A" / "G.ScandensA.stl" 

skull_name19 =  mesh_dir / "Full Skulls STL" / "C.ParvulusA.stl"
beak_name19 =  mesh_dir / "Cut Beaks B" / "C.ParvulusA_Cut.stl"
CH_name19 =  mesh_dir / "Convex Hulls A" / "C.ParvulusA_CH.stl"
WTR_name19 =  mesh_dir / "Watertight B" / "C.ParvulusA_WTR.stl" 
top_name19 =  mesh_dir / "Tops A" / "C.ParvulusA.stl" 
fit_name19 =  mesh_dir / "Fits A" / "C.ParvulusA.stl" 

skull_name20 =  mesh_dir / "Full Skulls STL" / "C.ParvulusC.stl"
beak_name20 =  mesh_dir / "Cut Beaks B" / "C.ParvulusC_Cut.stl"  
CH_name20 =  mesh_dir / "Convex Hulls A" / "C.ParvulusC_CH.stl"
WTR_name20 =  mesh_dir / "Watertight B" / "C.ParvulusC_WTR.stl" 
top_name20 =  mesh_dir / "Tops A" / "C.ParvulusC.stl" 
fit_name20 =  mesh_dir / "Fits A" / "C.ParvulusC.stl" 

skull_name21 =  mesh_dir / "Full Skulls STL" / "C2.FuscaA.stl"
beak_name21 =  mesh_dir / "Cut Beaks B" / "C2.FuscaA_Cut.stl"  
CH_name21 =  mesh_dir / "Convex Hulls A" / "C2.FuscaA_CH.stl"   
WTR_name21 =  mesh_dir / "Watertight B" / "C2.FuscaA_WTR.stl" 
top_name21 =  mesh_dir / "Tops A" / "C2.FuscaA.stl" 
fit_name21 =  mesh_dir / "Fits A" / "C2.FuscaA.stl" 

skull_name22 =  mesh_dir / "Full Skulls STL" / "C2.OlivaceaA.stl"
beak_name22 =  mesh_dir / "Cut Beaks B" / "C2.OlivaceaA_Cut.stl"  
CH_name22 =  mesh_dir / "Convex Hulls A" / "C2.OlivaceaA_CH.stl"  
WTR_name22 =  mesh_dir / "Watertight B" / "C2.OlivaceaA_WTR.stl"  
top_name22 =  mesh_dir / "Tops A" / "C2.OlivaceaA.stl" 
fit_name22 =  mesh_dir / "Fits A" / "C2.OlivaceaA.stl" 

skull_name23 =  mesh_dir / "Full Skulls STL" / "C2.OlivaceaB.stl"
beak_name23 =  mesh_dir / "Cut Beaks B" / "C2.OlivaceaB_Cut.stl"  
CH_name23 =  mesh_dir / "Convex Hulls A" / "C2.OlivaceaB_CH.stl"   
WTR_name23 =  mesh_dir / "Watertight B" / "C2.OlivaceaB_WTR.stl" 
top_name23 =  mesh_dir / "Tops A" / "C2.OlivaceaB.stl" 
fit_name23 =  mesh_dir / "Fits A" / "C2.OlivaceaB.stl" 

skull_name24 =  mesh_dir / "Full Skulls STL" / "C2.OlivaceaC.stl"
beak_name24 =  mesh_dir / "Cut Beaks B" / "C2.OlivaceaC_Cut.stl"  
CH_name24 =  mesh_dir / "Convex Hulls A" / "C2.OlivaceaC_CH.stl"   
WTR_name24 =  mesh_dir / "Watertight B" / "C2.OlivaceaC_WTR.stl" 
top_name24 =  mesh_dir / "Tops A" / "C2.OlivaceaC.stl" 
fit_name24 =  mesh_dir / "Fits A" / "C2.OlivaceaC.stl" 

skull_name25 =  mesh_dir / "Full Skulls STL" / "G.ConirostrisA.stl"
beak_name25 =  mesh_dir / "Cut Beaks B" / "G.ConirostrisA_Cut.stl"  
CH_name25 =  mesh_dir / "Convex Hulls A" / "G.ConirostrisA_CH.stl"   
WTR_name25 =  mesh_dir / "Watertight B" / "G.ConirostrisA_WTR.stl" 
top_name25 =  mesh_dir / "Tops A" / "G.ConirostrisA.stl" 
fit_name25 =  mesh_dir / "Fits A" / "G.ConirostrisA.stl" 

skull_name26 =  mesh_dir / "Full Skulls STL" / "G.ConirostrisB.stl"
beak_name26 =  mesh_dir / "Cut Beaks B" / "G.ConirostrisB_Cut.stl"  
CH_name26 =  mesh_dir / "Convex Hulls A" / "G.ConirostrisB_CH.stl"   
WTR_name26 =  mesh_dir / "Watertight B" / "G.ConirostrisB_WTR.stl" 
top_name26 =  mesh_dir / "Tops A" / "G.ConirostrisB.stl" 
fit_name26 =  mesh_dir / "Fits A" / "G.ConirostrisB.stl" 

skull_name27 =  mesh_dir / "Full Skulls STL" / "G.ConirostrisC.stl"
beak_name27 =  mesh_dir / "Cut Beaks B" / "G.ConirostrisC_Cut.stl"  
CH_name27 =  mesh_dir / "Convex Hulls A" / "G.ConirostrisC_CH.stl"  
WTR_name27 =  mesh_dir / "Watertight B" / "G.ConirostrisC_WTR.stl"  
top_name27 =  mesh_dir / "Tops A" / "G.ConirostrisC.stl" 
fit_name27 =  mesh_dir / "Fits A" / "G.ConirostrisC.stl" 

skull_name28 =  mesh_dir / "Full Skulls STL" / "G.SeptentrionalistA.stl"
beak_name28 =  mesh_dir / "Cut Beaks B" / "G.SeptentrionalistA_Cut.stl"  
CH_name28 =  mesh_dir / "Convex Hulls A" / "G.SeptentrionalistA_CH.stl"   
WTR_name28 =  mesh_dir / "Watertight B" / "G.SeptentrionalistA_WTR.stl" 
top_name28 =  mesh_dir / "Tops A" / "G.SeptentrionalistA.stl" 
fit_name28 =  mesh_dir / "Fits A" / "G.SeptentrionalistA.stl" 

skull_name29 =  mesh_dir / "Full Skulls STL" / "G.SeptentrionalistB.stl"
beak_name29 =  mesh_dir / "Cut Beaks B" / "G.SeptentrionalistB_Cut.stl"  
CH_name29 =  mesh_dir / "Convex Hulls A" / "G.SeptentrionalistB_CH.stl"   
WTR_name29 =  mesh_dir / "Watertight B" / "G.SeptentrionalistB_WTR.stl" 
top_name29 =  mesh_dir / "Tops A" / "G.SeptentrionalistB.stl" 
fit_name29 =  mesh_dir / "Fits A" / "G.SeptentrionalistB.stl" 

skull_name30 =  mesh_dir / "Full Skulls STL" / "P.CrassirostrisA.stl"
beak_name30 =  mesh_dir / "Cut Beaks B" / "P.CrassirostrisA_Cut.stl"  
CH_name30 =  mesh_dir / "Convex Hulls A" / "P.CrassirostrisA_CH.stl"   
WTR_name30 =  mesh_dir / "Watertight B" / "P.CrassirostrisA_WTR.stl" 
top_name30 =  mesh_dir / "Tops A" / "P.CrassirostrisA.stl" 
fit_name30 =  mesh_dir / "Fits A" / "P.CrassirostrisA.stl" 

skull_name31 =  mesh_dir / "Full Skulls STL" / "P.CrassirostrisB.stl"
beak_name31 =  mesh_dir / "Cut Beaks B" / "P.CrassirostrisB_Cut.stl"  
CH_name31 =  mesh_dir / "Convex Hulls A" / "P.CrassirostrisB_CH.stl"   
WTR_name31 =  mesh_dir / "Watertight B" / "P.CrassirostrisB_WTR.stl" 
top_name31 =  mesh_dir / "Tops A" / "P.CrassirostrisB.stl" 
fit_name31 =  mesh_dir / "Fits A" / "P.CrassirostrisB.stl" 

skull_name32 =  mesh_dir / "Full Skulls STL" / "P.CrassirostrisC.stl"
beak_name32 =  mesh_dir / "Cut Beaks B" / "P.CrassirostrisC_Cut.stl"  
CH_name32 =  mesh_dir / "Convex Hulls A" / "P.CrassirostrisC_CH.stl"  
WTR_name32 =  mesh_dir / "Watertight B" / "P.CrassirostrisC_WTR.stl"  
top_name32 =  mesh_dir / "Tops A" / "P.CrassirostrisC.stl" 
fit_name32 =  mesh_dir / "Fits A" / "P.CrassirostrisC.stl"  

skull_name33 =  mesh_dir / "Full Skulls STL" / "C2.FuscaB.stl"
beak_name33 =  mesh_dir / "Cut Beaks B" / "C2.FuscaB_Cut.stl"  
CH_name33 =  mesh_dir / "Convex Hulls A" / "C2.FuscaB_CH.stl" 
WTR_name33 =  mesh_dir / "Watertight B" / "C2.FuscaB_WTR.stl" 
top_name33 =  mesh_dir / "Tops A" / "C2.FuscaB.stl" 
fit_name33 =  mesh_dir / "Fits A" / "C2.FuscaB.stl" 

skull_name34 =  mesh_dir / "Full Skulls STL" / "P2.InornataA.stl"
beak_name34 =  mesh_dir / "Cut Beaks B" / "P2.InornataA_Cut.stl"  
CH_name34 =  mesh_dir / "Convex Hulls A" / "P2.InornataA_CH.stl" 
WTR_name34 =  mesh_dir / "Watertight B" / "P2.InornataA_WTR.stl" 
top_name34 =  mesh_dir / "Tops A" / "P2.InornataA.stl" 
fit_name34 =  mesh_dir / "Fits A" / "P2.InornataA.stl" 

beaks = np.array([beak_name1, beak_name2,beak_name3,beak_name4,beak_name5,beak_name6
         ,beak_name7, beak_name8,beak_name9,beak_name10,beak_name11,beak_name12,
         beak_name13, beak_name14,beak_name15,beak_name16,beak_name17,beak_name18
         ,beak_name19, beak_name20,beak_name21,beak_name22,beak_name23,beak_name24,
         beak_name25, beak_name26,beak_name27,beak_name28,beak_name29,beak_name30,
         beak_name31, beak_name32,beak_name33] )
    
hulls = np.array([CH_name1, CH_name2,CH_name3,CH_name4,CH_name5,CH_name6
         ,CH_name7, CH_name8,CH_name9,CH_name10,CH_name11,CH_name12,
         CH_name13, CH_name14,CH_name15,CH_name16,CH_name17,CH_name18
         ,CH_name19, CH_name20,CH_name21,CH_name22,CH_name23,CH_name24,
         CH_name25, CH_name26,CH_name27,CH_name28,CH_name29,CH_name30,
         CH_name31, CH_name32,CH_name33] )
    
watertights = np.array([WTR_name1, WTR_name2,WTR_name3,WTR_name4,WTR_name5,WTR_name6
         ,WTR_name7, WTR_name8, WTR_name9, WTR_name10, WTR_name11, WTR_name12,
         WTR_name13, WTR_name14, WTR_name15, WTR_name16, WTR_name17, WTR_name18
         , WTR_name19, WTR_name20, WTR_name21, WTR_name22, WTR_name23, WTR_name24,
         WTR_name25, WTR_name26, WTR_name27, WTR_name28, WTR_name29, WTR_name30,
         WTR_name31, WTR_name32, WTR_name33] )
    
top_names = np.array([top_name1, top_name2, top_name3, top_name4, top_name5,top_name6
         ,top_name7, top_name8, top_name9, top_name10, top_name11, top_name12,
         top_name13, top_name14, top_name15, top_name16, top_name17, top_name18
         , top_name19, top_name20, top_name21, top_name22, top_name23, top_name24,
         top_name25, top_name26, top_name27, top_name28, top_name29, top_name30,
         top_name31, top_name32, top_name33] )


fit_names = np.array([fit_name1, fit_name2,fit_name3,fit_name4,fit_name5,fit_name6
         ,fit_name7, fit_name8, fit_name9, fit_name10, fit_name11, fit_name12,
         fit_name13, fit_name14, fit_name15, fit_name16, fit_name17, fit_name18
         , fit_name19, fit_name20, fit_name21, fit_name22, fit_name23, fit_name24,
         fit_name25, fit_name26, fit_name27, fit_name28, fit_name29, fit_name30,
         fit_name31, fit_name32, fit_name33] )



    
    
ReducedN1 =  mesh_dir / "Watertight B" / "G.FortisA_WTR.stl"   
ReducedN2 =  mesh_dir / "Watertight B" / "G.MagnirostrisA_WTR.stl"   
ReducedN3 =  mesh_dir / "Watertight B" / "C.ParvulusB_WTR.stl"   
ReducedN4 =  mesh_dir / "Watertight B" / "C.PsittaculaA_WTR.stl"   
ReducedN5 =  mesh_dir / "Watertight B" / "G.FortisB_WTR.stl"   
ReducedN6 =  mesh_dir / "Watertight B" / "G.FuliginosaB_WTR.stl"   
ReducedN7 =  mesh_dir / "Watertight B" / "C.PallidusA_WTR.stl"   
ReducedN8 =  mesh_dir / "Watertight B" / "G.ScandensA_WTR.stl"   
ReducedN9 =  mesh_dir / "Watertight B" / "C.ParvulusA_WTR.stl"   

'''
face1 = np.loadtxt(mesh_dir / "Watertight B" / "G.FortisA_Cut_genus0_reduced_faces.txt", dtype=int) - 1
verts1 =  np.loadtxt(mesh_dir / "Watertight B" / "G.FortisA_Cut_genus0_reduced_vertices.txt")
face2 =  np.loadtxt(mesh_dir / "Watertight B" / "G.MagnirostrisA_Cut_genus0_reduced_faces.txt", dtype=int) - 1
verts2 =  np.loadtxt(mesh_dir / "Watertight B" / "G.MagnirostrisA_Cut_genus0_reduced_vertices.txt")
face3 = np.loadtxt(mesh_dir / "Watertight B" / "C.ParvulusB_Cut_genus0_reduced_faces.txt"   , dtype=int) - 1
verts3 = np.loadtxt( mesh_dir / "Watertight B" / "C.ParvulusB_Cut_genus0_reduced_vertices.txt" )
face4 =  np.loadtxt(mesh_dir / "Watertight B" / "C.PsittaculaA_Cut_genus0_reduced_faces.txt"   , dtype=int) - 1
verts4 =  np.loadtxt(mesh_dir / "Watertight B" / "C.PsittaculaA_Cut_genus0_reduced_vertices.txt"  )
face5 =  np.loadtxt( mesh_dir / "Watertight B" / "G.FortisB_Cut_genus0_reduced_faces.txt", dtype=int) - 1 
verts5 =  np.loadtxt( mesh_dir / "Watertight B" / "G.FortisB_Cut_genus0_reduced_vertices.txt")
face6 =  np.loadtxt(mesh_dir / "Watertight B" / "G.FuliginosaB_Cut_genus0_reduced_faces.txt"    , dtype=int) - 1
verts6 =  np.loadtxt(mesh_dir / "Watertight B" / "G.FuliginosaB_Cut_genus0_reduced_vertices.txt")
face7 =  np.loadtxt(mesh_dir / "Watertight B" / "C.PallidusA_Cut_genus0_reduced_faces.txt"   , dtype=int) - 1
verts7 =  np.loadtxt(mesh_dir / "Watertight B" / "C.PallidusA_Cut_genus0_reduced_vertices.txt"   )
face8 =  np.loadtxt(mesh_dir / "Watertight B" / "G.ScandensA_Cut_genus0_reduced_faces.txt"   , dtype=int) - 1
verts8 =  np.loadtxt(mesh_dir / "Watertight B" / "G.ScandensA_Cut_genus0_reduced_vertices.txt"  )
face9 =  np.loadtxt(mesh_dir / "Watertight B" / "C.ParvulusA_Cut_genus0_reduced_faces.txt"   , dtype=int) - 1
verts9 = np.loadtxt( mesh_dir / "Watertight B" / "C.ParvulusA_Cut_genus0_reduced_vertices.txt"  )
'''






'''
#=========================================================================================================================
# fits a parabola to the data
#========================================================================================================================= 
def fit_beak_shape(mesh_name):
    
    vertices, faces = get_connectivity(mesh_name) 
    
    x_pos = vertices[:,0]; y_pos = vertices[:,1]; z_pos = vertices[:,2]
    
    popt, pcov= op.curve_fit(beak_shape_form, np.vstack((y_pos, z_pos)), x_pos, p0=[1.5,0.4,1.5,0.4])
    
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    z = np.linspace(z_pos.min(), z_pos.max(), 100)
    y = np.linspace(y_pos.min(), y_pos.max(), 100)
    
   # z = np.linspace(0, 1, 100)
  #  y = np.linspace(0, 1, 100)
    Z = np.outer(z, np.ones_like(z))
    Y = np.outer(np.ones_like(y), y)
    X = beak_shape_form((Y, Z), *popt)
    
   
    
    ax.plot_surface(X, Y, Z)
    
    
    plt.show()
    #plt.plot(beak_shape_form((y_pos, z_pos), *popt), y_pos, z_pos, 'g--',
    #      label='fit: a=%5.3f, b=%5.3f, c=%5.3f, d=%5.3f' % tuple(popt))
    
    return popt, pcov
   
#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================    


#=========================================================================================================================
# functional form of beak shape
#========================================================================================================================= 
def beak_shape_form(data,a,b,c,d):
    y, z = data
    return -(((y-b)/a)**2+((z-d)/c)**2)+1.0
    


#=========================================================================================================================
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#=========================================================================================================================    

'''



















    
    
    
    
    
    
    